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ABSTRACT 

Supermassive stars (SMSs; > 10^ Mq) and their remnant black holes are promis- 
ing progenitors for supermassive black holes (SMBHs) observed in the early universe 
at z > 7. It has been postulated that SMSs forms through very rapid mass accretion 
onto a protostar at a high rate exceeding 0.01 Mq yr^^. According to recent studies, 
such rapidly accreting protostars evolve into "supergiant protostars", i.e. protostars 
consisting of a bloated envelope and a contracting core, similar to giant star. However, 
like massive stars as well as giant stars, both of which are known to be pulsationally 
unstable, supergiant protostars may also be also unstable to launch strong pulsation- 
driven outflows. If this is the case, the stellar growth via accretion will be hindered by 
the mass loss. We here study the pulsational stability of the supergiant protostars in 
the mass range < 10'^ Mq through the method of the linear perturbation analysis. 
We find that the supergiant protostars with il/* > 600 Mq and very high accretion rate 

Ma,cc ^1-0 Mq yr~^ are unstable due to the k mechanism. The pulsation is excited 
in the He"*" ionization layer in the envelope. Even under a conservative assumption 
that all the pulsation energy is converted into the kinetic energy of the outflows, the 
mass-loss rate is ^ 10~^ Mq yr"-'^, which is lower than the accretion rate by more 
than two orders of magnitude. We thus conclude that the supergiant protostars should 
grow stably via rapid accretion at least in the mass range we studied. As long as the 
rapid accretion is maintained in the later stage, protostars will become SMSs, which 
eventually produce seeds for the high-z SMBHs. 

Key words: stars: Population HI, protostars, oscillations, mass-loss - cosmology: 
theory - early Universe - galaxies: formation, nuclei 



1 INTRODUCTION 



reveal that 



super- 



Recent observations of high-z quasars 
massive black holes (SMBHs) of Mbh > lO'' Mq have 
already formed as early as the beginning of the universe 
< O.SGyr (e.g.. Fan 2006; Willott et al. 2007). A popular 
formation scenario of those SMBHs postulates that rem- 
nant stellar- mass BHs (Afgcod ~ 100 Mq) of Population 111 
(Pop HI) stars grow in mass via continuous mass accretion 
and merge (e.g., Haiman & Loeb 2001; Volonteri, Haardt & 
Madau 2003; Li et al. 2007). Given that the seed BHs grow 
at the Eddington mass-accretion rate MEdd ~ iEdd/ec^, 
where I/Edd is the Eddington luminosity, and e ~ 0.1 is 
the radiative efficiency, the growth time to lO'' Mq BHs 
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is ~ 0.051n(MBH/Ms<,cd) Gyr ~ 0.8 Gyr. Since this growth 
time is as long as the age of the universe at z ~ 7, where 
the most distant SMBH is observed (Mortlock et al. 2011), 
the stellar-mass seed is required to keep growing at least 
with the Eddington rate. However, recent studies show that 
this is unlikely as the accretion onto the BH, as well as onto 
the surrounding disk, is easily quenched by strong radiative 
feedback from the growing BH itself (Johnson & Bromm 
2007; Milosavljevic, Couch & Bromm 2009; Alvarez, Wise 
& Abel 2009; Park & Ricotti 2011; Park & Ricotti 2012; 
Tanaka, Perna & Haiman 2012). 

In an alternative scenario, formation of supermassive 
stars (SMSs; M* > lO'^ Mq) and their subsequent col- 
lapse directly to the BHs in the first galaxies {z > 10, 

> 10* K) has been envisaged (e.g., Bromm & Loeb 
2003; Begelman, Volonteri & Rees 2006; Lodato &i Natara- 
jan 2006). Here, primordial-gas clouds more massive than 
10^ Mq are supposed to contract monolithically to form 
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stars without strong fragmentation. Since rapid H2 cooling 
causes fragmentation of the primordial gas cloud, for the 
SMS formation, suppression of H2 formation is required by 
some means. Examples of such means are: the photodissoci- 
ation by far ultraviolet (FUV) radiation from nearby stars 
(e.g., Omukai 2001; Bromm & Loeb 2003; Omukai, Schnei- 
der & Haiman 2008; Regan & Haehnelt 2009a,b; Shang, 
Bryan, & Haiman 2010; Inayoshi & Omukai 2011; Agar- 
wal et a. 2012; Johnson, Dalla Vecchia & Khochfar 2012), 
and the collisional dissociation in dense and hot gas (In- 
ayoshi & Omukai 2012). The latter situation can be realized, 
for example, by the cold-accretion-flow shocks in the first 
galaxy formation. In both cases, the primordial gas collapses 
isothermally at T ~ 8000 K via H atomic cooling (Lya, two- 
photon, and H~ free-bound emission; Omukai 2001) and no 
major fragmentation is observed in numerical simulations 
during this phase (Bromm & Loeb 2003; Regan & Haehnelt 
2009a, b). 

This monolithic contraction of the cloud leads to a for- 
mation of a small protostar {r^ 0.01 Mq) at its center. 
The embryo protostar subsequently grows to a SMS via 
rapid accretion of the surrounding envelope. This process 
sounds similar to the case of ordinary Pop HI star for- 
mation, where H2 is the efficient coolant. However, there 
is an important difference; in the H atomic cooling case, 
the accretion rate onto the protostar is ~ 0.1 Mq yr~^, 
which is much higher than that in the ordinary Pop HI case 
(~ 10~^ Mq yr~^). This is due to the high temperature in 
the atomic cooling cloud (~ 8000 K) because the accretion 
rate is set by the temperature in the star-forming cloud as 
Mace ~ 10"^ Mq yr-^(r/600 K)^-^ (Shu 1977; Stabler et 
al. 1986). 

This rapid accretion with ~ 0.1 Mq yr^^ drasti- 
cally changes the protostellar evolution. Figure [T] shows 
the evolution of the radii of accreting protostars at dif- 
ferent accretion rates. In the ordinary Pop HI protostar 
case (Mace — 10~^ Mq yr~^), after the so-called adiabatic- 
accretion phase, where adiabatic heat input expands the star 
gradually with mass and the protostar starts to contract 
by losing its entropy via radiative diffusion (the Kelvin- 
HelmhoHz contraction) until the nuclear ignition occurs at 
the center. The protostar reaches the zero-age main sequence 
(ZAMS) stage at this point (Stahler, Palla & Salpeter 1986; 
Omukai & PaUa 2001, 2003). On the other hand, with the 
accretion rate as high as Maec ^ 0.1 Mq yr~^, the pro- 
tostar continues expanding without the KH contraction as 
recently found by Hosokawa, Omukai & Yorke (2012, here- 
after HOY12) (see Figure [l]). In such a star, while most of 
the interior material contracts, the outermost layer signifi- 
cantly swells up like a red-giant star (" supergiant protostar'^ 
phase). This is because the outer layer absorbs a part of 
the outward heat flux and obtains a very high specific en- 
tropy. Also in this case the contraction at the center ceases 
with the hydrogen ignition, but the envelope continuously 
expands with the increase of stellar mass. 

If rapid accretion at Aface ^ 0.1 Mq yr~^ is maintained, 
the stellar mass exceeds 10^ Mq within its lifetime. Such 
SMSs are general-relativistically unstable (e.g., Zel'dovich 
& Novikov 1971; Shapiro & Teukolsky 1983) and collapse 
as a whole to a BH (Shibata & Shapiro 2002), which can 
be a seed for the SMBHs residing in the early universe 
(2: > 7). With the stellar mass increasing, however, the stars 



become more radiation-pressure dominated and approach a 
marginally stable state. This may induce pulsational insta- 
bility of the massive stars and result in mass-loss from the 
surface. If this mass loss surpasses the accretion onto the 
star, the stellar growth will be terminated at that point. To 
see whether the SMS formations are indeed possible in spite 
of such mass loss, we examine the pulsational stability of 
the supergiant protostars in this paper. 

Pulsational stability of non-accreting Pop HI stars has 
been studied by Baraffe, Heger & Woosley (2001) and Sonoi 
& Umeda (2012) in the range 120 Mq ^ M. ^ 3 x 10^ Mq. 
They showed that those stars are unstable against pulsa- 
tion caused by the nuclear burning (the so-called e mech- 
anism), and that the resulting mass-loss rate is Mioss ^ 
10~^ Mq yr~^. Gamgami (2007) studied this mass-loss pro- 
cess from the Pop III stars using spherically symmetric hy- 
drodynamical simulations, and showed that pulsation accel- 
erates the surface material to the escape velocity and causes 
eruptive mass-loss for M* > 500 Mq . On the other hand, the 
Pop I red-giant stars are known to be pulsationally unstable 
by the opacity-driven mechanism (the so-called k mecha- 
nism, e.g., Li & Gong 1994; Heger et al. 1997), and the typ- 
ical mass-loss rate is ~ 10~^ Mq yr~^ (Yoon & Cantiello 
2010). With the Pop HI composition while having similar 
structure to the Pop I red-giants, the supergiant protostars 
can also be pulsationally unstable. 

In this paper, we study their stability by performing the 
linear stability analysis for the mass range M, < 10^ Mq, 
which has been calculated by HOY12. By estimating the 
mass-loss rate, we discuss whether supergiant protostars 
grow via accretion despite the pulsation-driven mass loss. 

The organization of this paper is as follows. In Section 
2, we introduce the method for the linear stability analysis 
against the pulsation and for the estimation of the mass- 
loss rates for unstable stars. In Section 3, we present our 
results and explain how the stability changes with the dif- 
ferent stellar masses and accretion rates. Finally, in Section 
4, we summarize our study and present our discussions. In 
Appendix, we describe the details (the basic equations and 
boundary conditions) of the linear stability analysis. 



2 STABILITY ANALYSIS 

We study the pulsational stability of protostars growing at 
constant accretion rates Maec = 10^'*, 0.03, 0.1, 0.3, and 
1.0 Mq yr~^, whose structures have been numerically cal- 
culated in our previous work (H0Y12). Figure[T]presents the 
evolution of the stellar radii for these rates. We apply the lin- 
ear stability analysis (see Appendix for the details) to stellar 
models either in the ZAMS (for Maec = 10"^Mq yr"^) or 
supergiant protostar (for higher accretion rates), indicated 
by the shaded areas in Figure [T] We consider the perturba- 
tions proportional to e"^', where cr = ctr + iai is the eigen 
frequency, ctr the frequency of the pulsation, and |ai| the 
growing or damping rate of the pulsation depending on its 
sign; the stars are stable (respectively unstable) if cri > O 
(respectively ai < 0). According to previous studies, mas- 
sive main-sequence Pop HI stars are unstable only under the 
radial perturbations (Baraffe et al. 2001; Sonoi & Umeda 
2012). We thus consider only the radial mode, hereafter, at 
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Figure 1. Evolution of the protostellar radius with various accre- 
tion rates Mace = 10-3, 0.03, 0.1, 0.3, and 1.0 Mq yr"! (taken 
from HOY12 with some modifications). In this paper, we analyze 
the pulsational stability of the stars located in the shaded zones; 
accreting ZAMS stars for Mace = lO""^ Mq yr" ^, and supergiant 
protostars for Mace > 10"^ Mq yr~l. 



which the supergiant protostars are also expected to be the 
most unstable. 

A useful quantity in the stability diagnosis is the work 
integral W (e.g.. Cox 1980; Unno et al. 1989), 

' ST* ( , d 



WiMr 



'^R Jo 



5R 



dM, 



dMr, (1) 



where Mr is the enclosed mass, T is the temperature, e is 
the nuclear energy generation rate per unit mass, Lrad is the 
radiative luminosity, and the symbols with S represent the 
Lagrange perturbations, where symbol K denotes the real 
part of the quantity indicated in the bracket. The work inte- 
gral has the physical meaning of the pulsation energy gained 
inside Mr in a single period. If the sign of the work integral 
is positive at the stellar surface, i.e., W{M,) > 0, the stars 
gain kinetic energy in each period and are unstable. The 
pulsation amplitude increases during the growth timescale 
of the instability crf^. If VK(M*) < 0, on the other hand, 
the pulsation damps inside the stars and are stable. The 
first term in the bracket on the right-hand side of equation 
(O, proportional to 5e, represents the driving of instability 
by the nuclear burning (i.e., the e mechanism). The second 
term is related to the radiative energy transport. In most 
cases, the radiative diffusion damps the pulsation, and thus 
the second term contributes to the stabilization. However, 
in the surface layer where the opacity changes remarkably, 
the energy flux transported via radiation can be absorbed 
and be converted into the pulsation energy by the k mech- 
anism. The growth (or damping) rate of the pulsation per 
single period ri = — cti/ctr is written as 

^ _ W{M,) 

(see Cox 1980; Unno et al. 1989), where 



(2) 



\^r\ dMr 



(3) 



is the pulsation energy, and ^r is the radial displacement of 
fluid elements from their equilibrium positions. 



In stars which are unstable under linear perturbations, 
the pulsation amplitude will grow to the non-linear regime. 
Such a strong pulsation is expected to cause a mass loss from 
the stellar surface (Appenzeller 1970a, b; Papaloizou 1973a, 
b). Appenzeller (1970a, b) studied the non-linear growth 
of the radial-pulsation instability for Pop I massive stars 
of M, = 130 and 270 Mq using one- dimensional hydro- 
dynamical calculations. He showed that after the pulsation 
enters the non-linear regime, the surface velocity reaches the 
speed of sound and weak shocks emerge just inside the pho- 
tosphere. The shocks recurrently propagate outward and ac- 
celerate the gas in the surface layer (e.g., Lamers & Cassinelli 
1999), generating mass shells exceeeding the escape velocity 
which is then lost from the star. 

In this paper, we evaluate the mass-loss rate following 
Baraffe et al. (2001) and Sonoi & Umeda (2012). As shown 
by Appenzeller (1970a, b), outflows are launched when the 
pulsation velocity at the surface reaches the speed of sound 
Cs. At this moment, the pulsation amplitude at the surface 
is 



CTR 



(4) 



Using this, we can estimate the pulsation energy Ew as well 
as the work integral W. Assuming that all the pulsation 
energy is converted into the kinetic energy of the outflows, 
the mass-loss rate Mioss can be obtained from the energy 
conservation: 



Mloss 



^WiM,) 



'2aiEw, 



(5) 



where Vcsc ~ {2GMt /Rt)^^^ is the escape velocity. 

Note that the assumption of energy conservation above 
is not always valid because some pulsation energy can be lost 
by radiative dissipation. In fact, Papaloizou (1973a, b) ob- 
tains a mass-loss rate lower than that of Appenzeller(1970a, 
b) by one order of magnitude by including this effect. The 
mass-loss rate derived below can thus be regarded as a con- 
servative upper limit. 



3 RESULTS 

In this Section, we describe the results for two different 
regimes of the accretion rate separately: (a) high accretion- 
rate cases (Afaec ^ 0.03 Mq yr~^), where the accreting stars 
become supergiant protostars, and (b) a low accretion-rate 
case (Mace = 10~3 Mq yr~^), where it reaches the ordinary 
ZAMS. The high-rate regime corresponds to the cases of the 
SMS formations, while the lower rates are expected in the 
ordinary Pop HI star formation. The latter results are pre- 
sented for comparison with the previous studies (Baraffe et 
al. 2001; Sonoi & Umeda 2012). Since we have found that 
accreting protostars are unstable only for the radial mode 
without nodes (fundamental mode or "F-mode"), we present 
the results for the F-mode below. 



3.1 High accretion-rate cases 

(Mace ^ 0.03 Mq yr~^): supergiant protostars 

We see here the high accretion-rate cases Mace ^ 
0.03 Mq yr~^, where the protostars grow in mass through 
the supergiant- protostar phase (Figure [TJ. We first explain 



© 0000 RAS, MNRAS 000, 000-000 



4 K. Inayoshi, T. Hosokawa, and K. Omukai 



0.6 



0.4 



0.2 













- 


\ T/Tc 








\p/f5, \ 
\ \ / 


/ Mr/M* 













10 



10 

radius r/R* 



10-' 



Figure 2. The interior structure of the accreting 10 Mq pro- 
tostar with Mace = 1.0 Mq yr"-*^ as a function of the relative 
radius r/ij*. The lines present the radial profiles of the enclosed 
mass (solid), density (dashed), and temperature (dotted), respec- 
tively. The enclosed mass is normalized by the stellar mass and 



others are normalized by their central values; pc 
and Tc = 2.1 X lO'^ K. The vertical line at r/i?,, : 
the inner boundary of the convective envelope. 



the case with the highest accretion rate A/a, 
and then the cases with the lower rates. 



= 0.2 g cm 
0.25 denotes 



1.0 M0 yr- 



3.1.1 Highest Accretion-Rate Case (Khcc = 1-0 Mq yr~^j 

Figure [5] shows the stellar interior structure when the stellar 
mass reaches 10^ Mq with Mace ~ 1.0 Mq yr~^. No convec- 
tive core develops in the interior since the hydrogen burning 
has not yet started. The star instead consists of a radiative 
core and an outer convective layer. Although the convective 
layer only constitutes 10% of the stellar mass, and the re- 
maining 90% is the radiative core, it covers a large portion of 
the radius. This structure consisting of the central core and 
bloated envelope is similar to that of red-giant stars. With 
Afacc = 1.0 Mq yr~^, the protostar reaches this structure at 
M* > 200 Mq. In this evolutionary stage, the stellar lumi- 
nosity is close to the Eddington value (L, ~ I/Edd oc M,), 
and the effective temperature remains almost constant at 
Teff ~ 5000 K due to the strong temperature-dependence 
of the H~ bound-free opacity. With these two conditions, 
the mass-radius relationship of the supergiant protostars is 
analytically written as 

1/2 



R, ~ 8.2 X lO"" R, 







103 Mq 



(6) 



which well agrees with the numerical results (H0Y12). 

Figure |3] shows the spatial distributions of the work in- 
tegral for the radial F-mode at the stellar masses of 300, 
500, and 10^ M!©- The work integral W changes remark- 
ably near the surface (< 3 x lO'' K), in particular, around 
4 X 10^ K, due to a opacity bump by the He"*" ionization. 
At 300 and 500 Mq, the work integrals are negative at the 
surface and the stars are stable. At M, = 10'^ on the 

other hand, the work integral at the stellar surface is pos- 
itive, i.e., the protostar is pulsationally unstable by the k 
mechanism excited in the He^ ionization layer. 



All the work integrals shown in Figure [3] are constant 
in the outer H" and He° ionization layers since the radiative 
energy transport is efficient enough there. The k mechanism 
neither excite nor damp the pulsation. This can be seen by 
comparing the following two timescales, the cooling time in 
the layer outside a radius r in the unperturbed state (ther- 
mal timescale; e.g., Sonoi & Shibahashi 2011) 



tth 



and the period of the pulsation 



tdyn — 



2tv 



(7) 



(8) 



The open circles in the Figure |3] indicate the transition 
points where the two timescales equal each other {tth ~ 
tdyn). Outside this point, tth is shorter than tdyn as the den- 
sity and the specific heat decrease outward. We call this 
region where tth < idyn as the non-adiabatic zone. The vari- 
ation of the work integral is almost zero there because the en- 
tropy is rapidly dissipated during a pulsation period. Thus, 
the surface value of the work integral, which determines the 
pulsational stability of the star, is fixed at the transition 
point to the non-adiabatic zone, where tth = tdyn- 

As seen in Figure O the surface value of the work in- 
tegral W{Mt) increases with the stellar mass and eventu- 
ally becomes positive for > 500 Mq: the protostar becomes 
pulsationally unstable. Since the work integral grows in the 
He"'' ionization layer inside the transition point but remains 
constant outside. This increase of the surface value W{Mt,) 
with the stellar mass can be understood by the concomitant 
outward-shift of the transition point, which in turn can be 
explained by comparing the two timescales; 



and 



tth oc 



f'dyn 



Rl 




(9) 



(10) 



near the surface. In deriving the equation ([9]), we used the 
fact that the term cpTp in equation ((T)) changes only slightly 
for T < 4 X 10^ K in the range 10^ Mq < M* < 10^ Mq. 
Note that the dynamical timescale (equation llOp has the 
same dependence as the free-fall timescale of the star. Elim- 
inating Rt and L» in equations (|9]) and (|10|l with equation 
((6)1 and L* ~ I/Edd, we obtain: 



td 



tth 



(11) 



the thermal timescale becomes longer with respect to the 
dynamical timescale near the surface with increasing stellar 
mass. In other words, the surface layer becomes more adia- 
batic: the non-adiabatic zone on the surface layer becomes 
thinner and the transition point moves outward as seen in 
Figure [3l As a result, the surface value of the work integral 
increases and the stars become more unstable as the stellar 
mass increase. 

The growth rate of the pulsation r] and the resulting 
mass-loss rate Mioss are shown in Figure |4] as a function 
of the stellar mass. At i\f, ~ 600 Mq, the star becomes 
pulsationally unstable and the mass loss rate increases with 
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Figure 3. Radial distributions of the work integral W (in an 
arbitrary unit) near the stellar surface (T < 3 X 10^ K) for Mace = 
1.0 Mq yr-l. The lines represent the M* = 300 (dotted), 500 
(dashed), and lO'^ Mq (solid) stars. The shaded zones denote the 
ionization layers of He+, He, and H from left to right. Open circles 
mark the transition points, where the thermal timescale is equal 
to the dynamical timescale, tn-^ = tdyn- 

the stellar mass thereafter. At Af, ~ lO'^ the mass- 

loss rate reaches 2 x 10"'^ Mq yr~^, two orders of magni- 
tude higher than that in the ZAMS case with accretion rate 
Mace = 10"^ Mq yr"^ (see Sec. 3.2 below). This is, how- 
ever, still lower than the accretion rate by a factor of 500. 
In the case with spherical symmetry, therefore, pulsation- 
driven outflow would be completely quenched by the rapid 
accretion. With some angular momentum, the accretion onto 
the star proceeds mostly through a circumstellar disk. In this 
case, the outflow escapes unhindered in the polar directions 
where the stellar surface is not covered by the accreting flow. 
We thus expect that the supergiant protostar loses some ma- 
terial via bipolar pulsation-driven outflows, while simultane- 
ously growing in mass due to a more rapid accretion from 
the disk. 

3.1.2 Variation with different accretion rates 

Next we see the cases with lower accretion rate 0.03 — 
0.3 Mq yr~^. Figure [5] presents the growth rate 77 for the 
radial F-mode in these cases as functions of the stellar mass. 
Roughly speaking, at a given stellar mass, the growth rate 
rj is higher for higher accretion rates. In our analysis, stars 
are unstable (i.e. ?7 > 0) only in the two highest accretion 
rate cases,: those with 1.0 Mq yr"^ for M* > 600 Mq and 
with 0.3 Mq yr"^ for M, > 900 Mq. 

This tendency of instability toward higher accretion 
rates can be understood again from the outward-shift of the 
transition point between the adiabatic and non-adiabatic 
zones inside the IIe+ ionization layer, which makes the sur- 
face value of the work integral M^(M*) higher (see Section 
I3.1.1[) . The behavior of work integral W is shown in Figure [6] 
for three accretion rates of 0.1, 0.3, and 1.0 Mq yr~^ at 
M* = 10'^ Mq. The ratio of the timescales tdyn/ith depends 
on the accretion rate Mace only through the stellar surface 
luminosity L* (see equations (6] [T] and [8] and note that the 
stellar radius is independent of Afacc). As shown in Figure [7] 



Figure 4. The growth/damping rate r]{= — cti/ctr) and the 
mass-loss rate Mioas as a function of stellar mass for Ajfacc = 
1.0 Mq yr~^. The left (right) vertical axis shows Mioss iri unit of 
10~^ Mq yr~^ (77 in unit of 10~^, respectively). In the shaded 
area the star is stable against the radial pulsation, i.e., r] < 0. 
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Figure 5. The growth/damping rate r] = — (Tj/itr as a func- 
tion of stellar mass for Mace = 1.0 (solid), 0.3 (long-dashed), 
0.1 (short-dashed), 0.03 Mq yr~^ (dash-dotted), respectively. In 
the shaded area the star is stable against the radial pulsation 
{rj < 0). The symbols on the lines show the models for which we 
analyze the stability. Large open triangles on the cases with 0.1 
and 0.03 Mq yr~^ indicate the onset of the hydrogen burning. 

(a) the surface luminosity L« and the ratio tdyn/^th is lower 
for higher Mace. In other words, the surface region becomes 
more adiabatic and the transition point moves closer to the 
surface at higher Aface, which makes the surface value of 
work integral higher as well. 

The above relation of L* and Aface can be understood in 
the following way. With lower Aface, the central part of the 
star has longer time to lose its entropy and the star takes 
a more centrally concentrated structure at a given stellar 
mass maintaining the same stellar radius (see Figure [7] b). 
As the radiative energy transport is efficient in the inner hot 
and dense part, such a star has larger radiative core. Since 
the luminosity grows proportionally to the enclosed mass 
Mr inside the radiative core but remains roughly constant 
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Figure 6. The same as Figure [S] but for the stellar models at 
10"^ Mq with three different accretion rates Mace = 1-0 Mq yr"-' 
(solid), 0.3 Mq yr-l (dashed), and 0.1 Mq yr"! (dotted). 

outside (see Figure [7| a), the large radiative core at low Mace 
results in high value of surface luminosity L,. 

Although the overall behavior of the growth rate rj 
shown in Figure [5] can be understood with the above con- 
siderations, 77 evolves in a somewhat complicated way at 
Afacc = 0.03 and 0.1 Mq yr~^; i.e., rj decreases with mass 
early in the evolution. The reason is as follows. As seen 
above, a supergiant protostar becomes more centrally con- 
centrated and thus more stable (i.e., lower rj) with the in- 
creasing mass. At the same time, however, there is also a 
destabilization effect with mass which is due to the shrink- 
ing of the non-adiabatic layer of the surface, as discussed in 
Sec. 13.1.11 These two effects compete each other. With the 
highest accretion rates of Mace = 0.3 and 1.0 Mq yr~^, the 
destabilization is more important, while in the cases with low 
accretion rate of Mace = 0.03 and 0.1 Mq yr~^, the stabi- 
lization first dominates until the onset of hydrogen burning, 
after which the central concentration remains almost con- 
stant and the stabilizing effect no longer operates. Thus, at 
this point r] begins to increase as seen in Figure (5] 



3.2 Lowest accretion-rate case 

(Afacc = 10"^ Mq yr~^ ): Accreting ZAMS stars 

Next, we will see the lowest accretion-rate case in our cal- 
culation with Mace ~ 10"'^ Mq yr~^, where the protostar 
reaches the ZAMS at M, — 50 Mq after the KH contraction 
(e.g., Omukai & Palla 2003). 

Figure [8] presents the interior structure of the protostar 
M« = 10^ M!q: the radial profiles of the mass, the temper- 
ature, the density, and the nuclear energy production rate. 
The central hydrogen burning via CN-cycle renders the 95% 
of the stellar mass (covering 60% in radius) to be convec- 
tive. The vertical line (dot-dashed) in Figure [8] indicates the 
boundary between the convective core and the outer radia- 
tive envelope. 

Figure [9] shows the radial distributions of the work in- 
tegral W, its derivative dW/dMr, its nuclear-energy pro- 
duction rate e, and its opacity k within this star. The e- 
mechanism drives the pulsational instability and thus the 
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Figure 7. Comparison of the interior structure of 10'^ Mq proto- 
stars with different accretion rates. The panels (a) and (b) present 
the radial distributions of the luminosity and enclosed mass, re- 
spectively. In the both panels, the red portions denote the con- 
vective zones. 
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Figure 8. The same as Figure 2, but for the lO'' Mq protostar 
with Mace = 10-3 Mq yr-i. The energy production rate due to 
the nuclear burning is also shown (dotted). The enclosed mass is 
normalized by the stellar mass, and others are by their central 
values: = 11 g cm-^, Te = 1.3 X 10** K, and Ce = 6.9 X 10^ 
erg g-^. The vertical line at r/ij* ~ 0.6 denotes the outer 
boundary of the convective core. 
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Figure 9. Radial distributions of several quantities within the 
Mt = 10^ M0 star with A'/acc = 10"^ Mq yr'^. The work 
integral W (solid line) and its derivative dW/dMr (long-dashed 
line) are presented in arbitrary units. The nuclear energy produc- 
tion rate e (normalized by the central value) and the opacity k 
(in cm^ g~^) are plotted with the short-dashed and dotted lines, 
respectively. 



Figure 10. The growth/damping rate r](= —ui/a^), and the 
mass loss rate Mioss as a function of stellar mass for Mace = 
10-3 M0 yr-l. The left (right) vertical axis shows Miogg in unit 
of 10~^ Mq yr" (r] in unit of IQ— ^, respectively). The protostar 
is stable {rj < 0) against the radial pulsation (F-mode) in the 
shaded area. Open triangles indicate the onset of the hydrogen 
burning. 



work integral W increases inside the convective core. On the 
other hand, the pulsation is slightly damped (i.e. dW/dMr < 
0) in the radiative layer because of the energy dissipation. 
The K mechanism does not work as the opacity is almost 
constant in the envelope due to high surface temperature 
lO^K). The small mass inside the surface region cannot 
totally damp the pulsation excited by the e mechanism. The 
star is thus unstable, i.e., VK(M*) > 0. 

Figure [10] presents the growth rate of the pulsation 77 
and the resulting mass-loss rate Mioss as a function of the 
stellar mass. After hydrogen ignition at M, — 50 Mq, the 
star remains stable until 140 Mq when the stabilization by 
radiative damping overcomes the pulsation by the e mech- 
anism. As the stellar mass increases, however, the star be- 
comes more radiation-pressure dominated and the average 
adiabatic exponent of the star Fi = {dlnp/dlnp)s ap- 
proaches the marginal gravitational stability value of 4/3. 
Consequently, the pulsation becomes increasingly strong in 
the central convective core and exceeds the radiative damp- 
ing effect (e.g.. Cox 1980; Shapiro & Teukolsky 1983). The 
star becomes unstable at 140 Mq and the growth rate of 
the pulsation increases thereafter. The mass-loss rate is typ- 
ically Mioss ^ - 10"^ Mq yr"^ (the dotted line in 
Figure [TO)) . Since this mass- loss rate is lower than the ac- 
cretion rate Mace = 10~3 Mq yr~^, the stellar growth via 
accretion would not be prevented by the pulsation-driven 
mass-loss as in the supergiant protostar cases. 

Baraffe et al. (2001) and Sonoi & Umeda (2012) also 
studied the pulsational instability of non-accreting mas- 
sive Pop III stars. For example, Sonoi & Umeda (2012) 
estimate the growth rate of the pulsation and mass-loss 
rate for a 500 Mq star as 77 = 3.02 x 10~* and Mioss = 
2.0x 10~^ Mq yr~ , respectively. Although their growth rate 
agrees well with our results, their mass-loss rate is higher 
than ours by a factor of four. This difference in mass-loss 
rates comes from the different values of the surface density 
between the accreting and non-accreting stars. With accre- 



tion, the surface density is higher than that without accre- 
tion. In this case, the gas pressure is relatively higher than 
the radiation pressure, i.e., higher /? = Pgas/ptot. This re- 
sults in a lower sound velocity Cs = \fVvpJp (oc /3~^''^) at 
the surface and a lower pulsation amplitude at the onset of 
the mass loss, which is proportional to the speed of sound 
(equation m . Therefore, the pulsation energy iJw and the 
mass-loss rate become lower in the accreting case than in 
the non-accreting case. 



4 CONCLUSION AND DISCUSSION 

In this paper, we have studied the pulsational stability 
of primordial protostars growing via very rapid accretion, 
(Mace ~ 0.1 Mq yr~^), through the method of the linear 
perturbation analysis, which is expected in the case of su- 
permassive star formation in the early universe. We have 
evaluated mass-loss rate if the protostar is pulsationally un- 
stable and examined whether the mass loss is strong enough 
to prevent the stellar growth via the accretion. We focused 
on early stellar evolution of M* < 10^ Mq, which has been 
studied in our recent work (HOY12). Our results are sum- 
marized as follows. 

First, we have studied the high accretion-rate cases with 
Mace 0.03 Mq yr~^, where the protostar has the a con- 
tracting core and a bloated envelope similar to a giant star 
[supergiant protostar; H0Y12). With low effective temper- 
ature Tcff ~ 5000 K, the supergiant protostar has the H 
and He ionization layers within its envelope. We have found 
that although pulsation is driven due to blocking of ra- 
diative flux at the opacity bump from the He"*" ionization 
(the so-called k mechanism), the supergiant protostars are 
pulsationally unstable only with the highest accretion rate 
~ 1.0 Mq yr~^ we studied. In the lower accretion-rate cases, 
the protostars are stable at least until M« ~ 10^ Even 
in the most unstable cases, the mass-loss rates are typically 
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^ 10"'^ Mq yr~^, which are lower than their accretion rates 
by more than two orders of magnitude. We thus conclude 
that the mass loss driven by pulsation does not prevent the 
growth of the supergiant protostar via rapid accretion. 

Next, for comparison with the previous studies, we 
have analyzed a lower accretion-rate case with Mace ~ 
10~^ Mq yr~^, which is expected in the ordinary Pop 
III star formation. In this case, the protostar reaches the 
ZAMS at M, ~ 50 Mq after the KH contraction stage 
(e.g., Omukai & Palla 2001, 2003). We have found that the 
protostars are unstable by the e mechanism in the range 
M« > 140 Mq, where a large part of the stellar interior 
is radiation-pressure dominated. Estimated mass-loss rate 
10~® — 10~^ Mq yr~^ is roughly consistent with the previ- 
ous results for non-accreting stars (Baraffe et al. 2001; Sonoi 
& Umeda 2012) although smaller by few factors because of 
the difference in the surface density due to the accretion. 

In this paper, we have limited our analysis to the mass 
range M, < 10'^ Mq due to the lack of stellar data in the 
higher mass range. We here speculate the later evolution 
based on the current results. Further studies on the proto- 
stellar evolution for M, > 10"^ Mq as well as on its pulsa- 
tional stability are thus awaited. If we linearly extrapolate 
the mass-loss rate in the case of Mace ~ 1-0 Mq yr~^ shown 
in Figure [4] to higher mass range, 

~ 5.0 X 10-* (^Y^j^ - Mq yr-^ (12) 

the mass loss catches up with the accretion at M* ~ 
2 X 10^ M!©- At this point, the growth of the protostar via 
accretion possibly halts and the final mass is set. However, 
because of our assumption that all the pulsation energy is 
converted to the kinetic energy of the outflows, the mass-loss 
rate by equation (|12|) should be regarded as an upper limit 
(e.g., Papaloizou 1973a, b). Furthermore, the mass loss of 
the Pop I red-giant stars are usually driven by the radiation 
pressure exerted on dust grains formed in the cool envelope 
(e.g., Willson 2000). Without the dust as in our case, accel- 
eration of the outflows could be more inefhcient. We expect 
that, with such rapid accretion, the final stellar mass can 
exceed ~ 10^ Mq despite the pulsation-driven mass loss. 

Although mass-loss rate is much lower than the accre- 
tion rate until M« = 10^ Mq as studied in this paper, these 
values could be comparable to A/* > 10^ 1^0 ■ Since the ac- 
cretion of gas with some angular momentum onto the star 
proceeds via a circumstellar disk, the outflows would escape 
most easily in the polar regions, where the density is rela- 
tively low. The dynamical interaction between the inflows 
and outflows needs to be studied in detail to determine the 
exact value of the stellar final mass. 

Our estimate of the mass-loss rate is based on the previ- 
ous works (e.g., Appenzeller 1970a, b; Papaloizou 1973a, b), 
in which the non-linear development of pulsation for non- 
accreting main-sequence stars is studied numerically. For 
the accreting stars, however, we have a very limited knowl- 
edge on the non-linear behavior of pulsation (e.g., Gamgami 
2007). More detailed studies on this issue by radiative hy- 
drodynamical simulations is awaited. 

So far, we have only considered stars forming from the 
metal-free {Z = 0) gas. However, SMSs could potentially be 
formed from the gas slightly polluted with heavy elements, if 
that is below the critical amount Zcr, i.e., ^ Zq with- 



out dust grains, or ~ 10~^ Zq with dust grains (Omukai 
et al. 2008; Inayoshi & Omukai 2012). Metal enrichment 
lowers the central temperature of a star by enhancing the 
energy production efficiency by nuclear fusion and also cre- 
ates another opacity bump near the surface which make the 
pulsation stronger via the e and k mechanisms, respectively. 
However, for the n mechanism, for which the supergiant pro- 
tostars are unstable, this effects becomes important only 
with metallicity higher than 2 x 10"'^ Zq in the case of 
non-accreting stars (Baraffe et al. 2001), which is higher 
than the critical value Za- We thus speculate that even if 
small amount of metals below Za are present, the pulsa- 
tional stability of supergiant protostars should be similar to 
the zero-metallicity case studied above. 

Finally, we discuss the validity of the frozen-in approx- 
imation of convective energy flux used in our analysis (see 
Appendix A), where perturbations of the convective flux is 
neglected. At the present time, this is a widely-used approx- 
imation due to our limited knowledge on the interaction 
between the convective and pulsational motions. Although 
some other models including this effect have been proposed 
(e.g., Unno 1967; Gough 1977; Unno et al. 1989; Dupret et 
al. 2005), they rely on the still-developing time-dependent 
convection theories, which require different assumptions 
depending on modeling, beyond the classical mixing-length 
theory (Bohm-Vitense 1958). Recent results by Penev, 
Barranco & Sasselov (2009) and Shiode, Quataert & Arras 
(2012), who studied this interaction numerically, showed 
that the convective damping weakens the pulsation by 
the e mechanism, but does not influence through the k 
mechanism. Therefore, protostars with modest accretion 
rate Mace — 10~^ Mq yr~^, which are unstable by the 
e mechanism (Sec. 13. 2|) . can be somewhat stabilized by 
this convective damping. On the other hand, we speculate 
that this would not signiflcantly affect the pulsation in 
supergiant protostars, which is driven by the k mechanism. 
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APPENDIX A: LINEAR PERTURBATION 
ANALYSIS METHOD 

In this appendix, we describe our method of the linear per- 
turbation analysis of the stellar pulsation stability (e.g.. Cox 
1980; Unno et al. 1989). The basic equations governing the 
stellar structure are 

^ + ipv) = 0, (Al) 
^ + (v- V)v = -ivp- V$, (A2) 
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AnGp, 



dS_ 

dt 



+ (v- V)^ 



-V • F, 



3kp 



(A3) 



(A4) 



(A5) 



where p is the density, v the velocity, p the pressure, $ the 
gravitational potential, T the temperature, S the specific en- 
tropy, e the nuclear-energy generation rate per unit mass, k 
the opacity, F the total energy flux, which is the sum of the 
radiative flux Frad and convective flux Fcoiw, G the gravi- 
tational constant, c the speed of light, and a the radiation 
constant. The radial mode, i.e., perturbations which radially 
oscillate with an eigen frequency a, is studied. 

We consider the radial displacement of fluid elements 
from the equilibrium positions in the form £^r{T,t) = 
^r(r)e"^*. We define the resulting Euler perturbation of a 
physical quantity Q as Q' = Q{r,t) — Qo{r,t), where Qo is 
the value in the unperturbed state. We also use the Lagrange 
perturbation 5Q = Q{r + ^r, i) — Qo{r,t) for some physical 
quantities. The linearized equations (|A1|) - (|A5I) with the per- 



turbations Q'{r,t) 
written as 



?'(r)e""' and 5Q(r,t) = 5Q(r)e" 



^ d, 2c ^ 9 c , P' 

Cr-J - H J = -"T — , 

ar c| pc| cp 



(A6) 



+ ^p' + {N'^a')& + ^=gvT-, (A7) 
— ^ ar Cp 



1 rfp^ _l_ _9__^, , ,^r2 _2 
p dr pc? 



this relation, we obtain four linear ordinary first-order dif- 
ferential equations for four variables ^r, p' , SS, and 5Lrad. 
Here, we impose the following boundary conditions: 



dr \ r J ^' dr 



5Lr 



= (r = 0), (A12) 



dr 



(r = 



(A13) 



from the regularity of the perturbations at the center and 
surface, and 

^-4f {r^R,), (A14) 

which guarantees outward propagation of the energy flux at 
the surface (e.g.. Cox 1980; Saio, Winget & Robinson 1983). 
In this system of the differential equations and boundary 
conditions, the normalization of the variables ^r, p' , SS, 
SLrad still remains as a degree of freedom. We solve the 
system as an eigenvalue problem by formally introducing a 
differential equation for the eigenvalue a, 

The whole system here is the five first-order differential 
equations for £,r,p , SS, 5I/rad, and a with the four boundary 
conditions and one normalization condition. We set the arbi- 
trary normalization condition at the surface, (,r{r = R*) — 
R*. We obtain numerical solutions of the eigen functions 
and eigenvalue using the relaxation method (e.g., Unno et 
al. 1989). 



1 d / 2rf*'\ . ^ I P 1^" ^ \ . ^ <>S 

— -r- r^-r- - 47rGp + = -47rGp«T — , 

r-' dr \ dr J \ pci q J cp 

(A8) 



iaTSS = Se 



\pct g 

dSLmd 



dMr 



5L 



Sk ST ^ 
7 = +4^+4^ 



rf(4f)/rflnr 
dluT/dlnr- ' 



(A9) 



(AlO) 



where Ca (= vfTp/p) is the sound velocity, Fi = 
(91np/c'lnp)s the adiabatic exponent, r the radius, g the 
gravitational acceleration, Lrad the radiative luminosity, Mr 
the enclosed mass, cp — T{dS/dT)p the isobaric specific 
heat, vt = — (91np/91nr)p, and A^^ = ~'g{dlnp/dr + g/c^) 
the Brunt- Vaisara frequency. In the above equations, for 
simplicity, a physical quantity "Q" indicates its value in the 
unperturbed state instead of Qo- 

Note that, in equations (|A9|) and (|A10|) . we ignore the 
perturbation of the convective energy flux, i.e., 5Fconv ~ 0. 
This so-called "frozen-in" approximation has been widely 
used in analyzing the pulsational stability of stars (Baraffe 
et al. 2001 and Sonoi & Umeda 2012). To facilitate the com- 
parison between their results, we also adopt this assumption 
here (see Section |4] for more discussions). 

From equations HA6[) and HA8[I and the regularity of 
at the center. 



d$' 
dr 



+ ATxGp^r = 0. 



(All) 



Ehminating the term d^' /dr in equations (|A7I) - (|A10|I with 
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